#
# ANES2000.r
#
rm(list=ls(all=TRUE))
#
library(haven)
library(plyr)
library(dplyr)
library(questionr)
library(statar)
library(ggplot2)
library(viridis)
library(rjags)
library(apsrtable)
library(mokken)
library(readstata13)
library(basicspace)
library(gridExtra)
library(psych)
library(interflex)
library(ggpubr)
#
setwd("c:/")
data <- read_dta("Dropbox/Networks/data/anes2000TS.dta")
#
attach(data,warn.conflicts = FALSE)
#
rescale01 <- function(x){(x-min(x, na.rm=TRUE))/(max(x, na.rm=TRUE)-min(x, na.rm=TRUE))}
higher.conservative <- function(x){
	if (cor(x, ideology, use="complete") < 0) {
		(x <- -1 * x + max(x, na.rm=TRUE))}
	else{x <- x}}
#
#  PARTY ID (1 = Democrat, 2 = Independent, 3 = Republican)
#
partyid <- NA
partyid[V000523==0 | V000523==1 | V000523==2] <- 1
partyid[V000523==3] <- 2
partyid[V000523==4 | V000523==5 | V000523==6] <- 3
#
partyid7 <- V000523 + 1
partyid7[partyid7 < 1 | partyid7 > 7] <- NA
#
partisanextremity <- abs(partyid7 - 4)
#
#
# PRESIDENTIAL VOTE
#
presvote <- NA
presvote[V001249==1] <- "Gore"
presvote[V001249==3] <- "Bush"
presvote[V001249==2 | V001249==4 | V001249==5 | V001249==6 | V001249==7| V001249==8] <- "Other/not sure"
presvote <- as.factor(presvote)
#
#
# DEMOGRAPHICS
#
education <- V000913
education[education < 1 | education > 7] <- NA
education <- rescale01(education)
#
income <- V000997
income <- rescale01(income)
#
age <- V000908
age <- rescale01(age)
#
female <- NA
female[V001029==1] <- 0
female[V001029==2] <- 1
#
black <- NA
black[V001006a!=10] <- 0
black[V001006a==10] <- 1
#
latino <- NA
latino[V001006a!=40] <- 0
latino[V001006a==40] <- 1
#
unionmember <- NA
unionmember[V000990==5] <- 1
unionmember[V000990==1] <- 2
#
bornagain <- NA
bornagain[V000903==5] <- 1
bornagain[V000903==1] <- 2
bornagain <- (bornagain-min(bornagain, na.rm=TRUE))/(max(bornagain, na.rm=TRUE)-min(bornagain, na.rm=TRUE))
#
importancereligion <- NA
importancereligion[V000872==5] <- 1
importancereligion[V000872==1 & V000873==1] <- 2
importancereligion[V000872==1 & V000873==3] <- 3
importancereligion[V000872==1 & V000873==5] <- 4
importancereligion <- (importancereligion-min(importancereligion, na.rm=TRUE))/(max(importancereligion, na.rm=TRUE)-min(importancereligion, na.rm=TRUE))
#
churchattendance <- NA
churchattendance[V000877==5] <- 1
churchattendance[V000877==1 & V000879==4] <- 2
churchattendance[V000877==1 & V000879==3] <- 3
churchattendance[V000877==1 & V000879==2] <- 4
churchattendance[V000877==1 & V000879==1 & V000880==1] <- 5
churchattendance[V000877==1 & V000879==1 & V000880==2] <- 6
churchattendance <- (churchattendance-min(churchattendance, na.rm=TRUE))/(max(churchattendance, na.rm=TRUE)-min(churchattendance, na.rm=TRUE))
#
frequencyprayer <- NA
frequencyprayer[V000874==5] <- 1
frequencyprayer[V000874==4] <- 2
frequencyprayer[V000874==3] <- 3
frequencyprayer[V000874==2] <- 4
frequencyprayer[V000874==1] <- 5
frequencyprayer <- (frequencyprayer-min(frequencyprayer, na.rm=TRUE))/(max(frequencyprayer, na.rm=TRUE)-min(frequencyprayer, na.rm=TRUE))
#
religiositymat <- cbind(bornagain, importancereligion, churchattendance, frequencyprayer)
psych::alpha(religiositymat)
religiosity.scores <- rowMeans(religiositymat, na.rm=TRUE)
religiosity.binary <- ntile(religiosity.scores, 3)
#
#rr1 <- CC16_422c
#rr1[rr1 < 1 | rr1 > 5] <- NA
#
#rr2 <- CC16_422d
#rr2[rr2 < 1 | rr2 > 5] <- NA
#
#rr3 <- CC16_422e
#rr3[rr3 < 1 | rr3 > 5] <- NA
#rr3 <- -1 * (rr3 - 6)
#
#rr4 <- CC16_422f
#rr4[rr4 < 1 | rr4 > 5] <- NA
#rr4 <- -1 * (rr4 - 6)
#
#rrmat <- cbind(rr1, rr2, rr4)
#psych::alpha(rrmat)
#racialresentment.scores <- rowMeans(rrmat, na.rm=TRUE)
#racialresentment.binary <- ntile(racialresentment.scores, 3)
#
#
# COUNTY-LEVEL INFORMATION
#
countydata <- read_dta("Dropbox/Networks/data/anes2000_county_update.dta")
#
#county.ideology <- mrp_ideology_mean
#
county.gop2000 <- countydata$prop_gop_2000
#
#county.delta.gop2012gop2016 <- per_gop_2016 - per_gop_2012
#
#county.trumpprimary2016 <- trump_2016primary
#
#county.sandersprimary2016 <- sanders_2016primary
#
#
# SUMMATED IDEOLOGICAL SCALE
#
ideology <- NA
ideology[V000446==1] <- 1
ideology[V000446==2] <- 2
ideology[V000446==3] <- 3
ideology[V000446==4] <- 3
ideology[V000446==5] <- 3
ideology[V000446==6] <- 4
ideology[V000446==7] <- 5
#
govtservices <- V000545
govtservices[govtservices < 1] <- NA
govtservices[govtservices > 7] <- NA
#
govtservices.branch <- V000549
govtservices.branch[govtservices.branch < 1] <- NA
govtservices.branch[govtservices.branch > 5] <- NA
#
defensespending <- V000581
defensespending[defensespending < 1] <- NA
defensespending[defensespending > 7] <- NA
#
defensespending.branch <- V000586
defensespending.branch[defensespending.branch < 1] <- NA
defensespending.branch[defensespending.branch > 5] <- NA
#
healthinsurance <- V000609
healthinsurance[healthinsurance < 1] <- NA
healthinsurance[healthinsurance > 7] <- NA
#
healthinsurance.branch <- V000613
healthinsurance.branch[healthinsurance.branch < 1] <- NA
healthinsurance.branch[healthinsurance.branch > 5] <- NA
#
guarjobs <- V000615
guarjobs[guarjobs < 1] <- NA
guarjobs[guarjobs > 7] <- NA
#
guarjobs.branch <- V000619
guarjobs.branch[guarjobs.branch < 1] <- NA
guarjobs.branch[guarjobs.branch > 5] <- NA
#
aidtoblacks <- V000641
aidtoblacks[aidtoblacks < 1] <- NA
aidtoblacks[aidtoblacks > 7] <- NA
#
aidtoblacks.branch <- V000644
aidtoblacks.branch[aidtoblacks.branch < 1] <- NA
aidtoblacks.branch[aidtoblacks.branch > 7] <- NA
#
environmentjobs <- V000708
environmentjobs[environmentjobs < 1] <- NA
environmentjobs[environmentjobs > 7] <- NA
#
environmentjobs.branch <- V000712
environmentjobs.branch[environmentjobs.branch < 1] <- NA
environmentjobs.branch[environmentjobs.branch > 5] <- NA
#
womensrole.scale <- V000755
womensrole.scale[womensrole.scale < 1] <- NA
womensrole.scale[womensrole.scale > 7] <- NA
#
womensrole.branch <- V000759
womensrole.branch[womensrole.branch < 1] <- NA
womensrole.branch[womensrole.branch > 5] <- NA
#
abortion.scale <- V000694
abortion.scale[abortion.scale < 1] <- NA
abortion.scale[abortion.scale > 4] <- NA
#
gaymilitary <- V000727
gaymilitary[gaymilitary==4] <- 3
gaymilitary[gaymilitary==5] <- 4
gaymilitary[gaymilitary < 1] <- NA
gaymilitary[gaymilitary > 4] <- NA
#
immigrationlevel <- V000510
immigrationlevel[immigrationlevel < 1] <- NA
immigrationlevel[immigrationlevel > 5] <- NA
#
#
issues <- cbind(govtservices, govtservices.branch,
	defensespending, defensespending.branch,
	healthinsurance, healthinsurance.branch,
	guarjobs, guarjobs.branch,
	aidtoblacks, aidtoblacks.branch,
	environmentjobs, environmentjobs.branch,
	womensrole.scale, womensrole.branch,
	abortion.scale,
	gaymilitary,
	immigrationlevel)
#
issues <- apply(issues, 2, rescale01)
issues <- apply(issues, 2, higher.conservative)
#
policypreferences <- scale(rowMeans(issues, na.rm=TRUE))
#
#
#
#
# NUMBER OF DISCUSSANTS
#
#numberdiscussants <-  NA
#numberdiscussants[legitimate_name_1==0] <- 0
#numberdiscussants[!is.na(UCD401_1) & legitimate_name_1==1] <- 1
#numberdiscussants[!is.na(UCD401_2) & legitimate_name_1==1] <- 2
#numberdiscussants[!is.na(UCD401_3) & legitimate_name_1==1] <- 3
#numberdiscussants[is.na(numberdiscussants)] <- 0
#
#
#  DISCUSSANT VOTE (1 = Clinton, 2 = Trump, 3 = Other/Not Sure, 0 = Missing)
#
#discussant.1.vote <- NA
#discussant.1.vote[UCD405A==1 & numberdiscussants >= 1] <- 1
#discussant.1.vote[UCD405A==2 & numberdiscussants >= 1] <- 2
#discussant.1.vote[(UCD405A==3 | UCD405A==5)  & numberdiscussants >= 1] <- 3
#discussant.1.vote[is.na(discussant.1.vote)] <- 0
#
#discussant.2.vote <- NA
#discussant.2.vote[UCD405B==1 & numberdiscussants >= 2] <- 1
#discussant.2.vote[UCD405B==2 & numberdiscussants >= 2] <- 2
#discussant.2.vote[(UCD405B==3 | UCD405B==5) & numberdiscussants >= 2] <- 3
#discussant.2.vote[is.na(discussant.2.vote)] <- 0
#
#discussant.3.vote <- NA
#discussant.3.vote[UCD405C==1  & numberdiscussants >= 3] <- 1
#discussant.3.vote[UCD405C==2  & numberdiscussants >= 3] <- 2
#discussant.3.vote[(UCD405C==3 | UCD405C==5)  & numberdiscussants >= 3] <- 3
#discussant.3.vote[is.na(discussant.3.vote)] <- 0
#
#number.clinton.voters <- as.numeric(discussant.1.vote==1) + as.numeric(discussant.2.vote==1) + as.numeric(discussant.3.vote==1)
#number.trump.voters <- as.numeric(discussant.1.vote==2) + as.numeric(discussant.2.vote==2) + as.numeric(discussant.3.vote==2)
#number.other.voters <- as.numeric(discussant.1.vote==3) + as.numeric(discussant.2.vote==3) + as.numeric(discussant.3.vote==3)
#
prop.gore.voters <- countydata$dgoreprop
prop.bush.voters <- countydata$dbushprop
#
wtd.table(presvote, prop.gore.voters, weights=countydata$weight)
wtd.table(presvote, prop.bush.voters, weights=countydata$weight)
#
#
# POLITICAL SOPHISTICATION
#
voted1996 <- NA
voted1996[V000303==5] <- 0
voted1996[V000303==1] <- 1
#
voted2000 <- NA
voted2000[V001241==1] <- 0
voted2000[V001241==2] <- 0
voted2000[V001241==3] <- 0
voted2000[V001241==4] <- 1
#
campaigninterest <- NA
campaigninterest[V000301==5] <- 0
campaigninterest[V000301==3] <- 0.5
campaigninterest[V000301==1] <- 1
#
knowledge.trentlott <- NA
knowledge.trentlott[V001447==5] <- 0
knowledge.trentlott[V001447==8] <- 0
knowledge.trentlott[V001447==1] <- 1
#
knowledge.rehnquist <- NA
knowledge.rehnquist[V001450==5] <- 0
knowledge.rehnquist[V001450==8] <- 0
knowledge.rehnquist[V001450==1] <- 1
#
knowledge.blair <- NA
knowledge.blair[V001453==5] <- 0
knowledge.blair[V001453==8] <- 0
knowledge.blair[V001453==1] <- 1
#
knowledge.reno <- NA
knowledge.reno[V001456==5] <- 0
knowledge.reno[V001456==8] <- 0
knowledge.reno[V001456==1] <- 1
#
knowledge.controlhousesenate <- NA
knowledge.controlhousesenate[V001356==1] <- 0
knowledge.controlhousesenate[V001356==8] <- 0
knowledge.controlhousesenate[V001357==1] <- 0
knowledge.controlhousesenate[V001357==8] <- 0
knowledge.controlhousesenate[V001356==5 & V001357==5] <- 1
#
knowledge.controlhouse <- NA
knowledge.controlhouse[V001356==1] <- 0
knowledge.controlhouse[V001356==5] <- 1
knowledge.controlhouse[V001356==8] <- 0
#
knowledge.controlsenate <- NA
knowledge.controlsenate[V001357==1] <- 0
knowledge.controlsenate[V001357==5] <- 1
knowledge.controlsenate[V001357==8] <- 0
#
demparty <- V001382
demparty[demparty < 1 | demparty > 7] <- NA
#
repparty <- V001383
repparty[repparty < 1 | repparty > 7] <- NA
#
knowledge.ideologyparties <- NA
knowledge.ideologyparties[V001382==8 | V001382==9] <- 0
knowledge.ideologyparties[V001383==8 | V001383==9] <- 0
knowledge.ideologyparties[demparty >= repparty] <- 0
knowledge.ideologyparties[demparty < repparty] <- 1
#
politicalsophisticationmat <- cbind(knowledge.controlhouse,
	knowledge.controlsenate,
	knowledge.ideologyparties)
#
psych::alpha(politicalsophisticationmat)
politicalsophistication.scores <- rowMeans(politicalsophisticationmat, na.rm=TRUE)
politicalsophistication.binary <- NA
politicalsophistication.binary[politicalsophistication.scores < 1] <- "Low"
politicalsophistication.binary[politicalsophistication.scores >= 1] <- "High"
politicalsophistication.binary <- factor(politicalsophistication.binary, levels=c("Low","High"))
#
table(rowSums(is.na(politicalsophisticationmat)))
#
#library(foreign)
#write.dta(data.frame(politicalknowledgescores=politicalsophistication.scores), file="c:/Dropbox/Networks/analysis/ANES2000knowledge.dta")
#
networkmix <- NA
networkmix[prop.gore.voters==1] <- "Only Gore"
networkmix[prop.gore.voters > 0 & prop.gore.voters < 1 & prop.bush.voters==0] <- "Gore + not sure"
networkmix[prop.gore.voters > 0 & prop.bush.voters > 0] <- "Gore + Bush"
networkmix[prop.gore.voters==0 & prop.bush.voters > 0 & prop.bush.voters < 1] <- "Bush + not sure"
networkmix[prop.bush.voters==1] <- "Only Bush"
#
networkmix <- factor(networkmix, levels=c("Only Gore", "Gore + not sure", "Gore + Bush", "Bush + not sure", "Only Bush"))
#
#
round((table(networkmix) / sum(table(networkmix))), 2)
table(networkmix, politicalsophistication.binary)
colSums(table(networkmix, politicalsophistication.binary))
round((table(networkmix, politicalsophistication.binary) /
	colSums(table(networkmix, politicalsophistication.binary))), 2)
#
round((wtd.table(networkmix, weights=countydata$weight) / sum(wtd.table(networkmix, weights=countydata$weight))), 2)
wtd.table(networkmix, politicalsophistication.binary, weights=countydata$weight)
colSums(wtd.table(networkmix, politicalsophistication.binary, weights=countydata$weight))
round((wtd.table(networkmix, politicalsophistication.binary, weights=countydata$weight) /
	colSums(wtd.table(networkmix, politicalsophistication.binary, weights=countydata$weight))), 2)
#
#
countytype <- NA
countytype[county.gop2000 < .4] <- "Blue"
countytype[county.gop2000 >= .4 & county.gop2000 <= .6] <- "Purple"
countytype[county.gop2000 > .6] <- "Red"
#
weight <- countydata$weight
df <- data.frame(presvote, networkmix, politicalsophistication.binary, weight, countytype)
df2 <- na.omit(df)
#
#df2$primary2016vote <- factor(primary2016vote[complete.cases(df)], levels=c("Clinton", "Sanders", "Other GOP", "Trump"))
#
df3 <- ddply(df2,.(presvote),
    function(x) with(x,
      data.frame(100*wtd.table(networkmix, weights=weight)/sum(wtd.table(networkmix, weights=weight)))))
#
names(df3)[names(df3)=="Var1"] <- "Network"
names(df3)[names(df3)=="Freq"] <- "Percentage"
#
#
pdf("Dropbox/networks/images/anes2000_fig1.pdf", height=7, width=7)
#
ggplot(df3, aes(x=presvote, y=Percentage, fill=Network)) + ylim(0,100.1) +
geom_bar(stat="identity") + scale_fill_viridis(discrete=TRUE) + xlab("2016 presidential vote intention") +
ggtitle("Composition of voters' discussion networks\n2016 Cooperative Congressional Election Study (weighted)") +
theme(plot.title = element_text(hjust = 0.5))
#
dev.off()
#
#
#
# Political sophistication
#
df5 <- ddply(df2,.(presvote, politicalsophistication.binary),
    function(x) with(x,
      data.frame(100*wtd.table(networkmix, weights=weight)/sum(wtd.table(networkmix, weights=weight)))))
#
names(df5)[names(df5)=="politicalsophistication.binary"] <- "sophistication"
names(df5)[names(df5)=="Var1"] <- "Network"
names(df5)[names(df5)=="Freq"] <- "Percentage"
#
df5 <- na.omit(df5)
#df5$sophistication <- factor(df5$sophistication, levels=c("Low", "Mid", "High"))
#
#
pdf("Dropbox/networks/images/anes2000_fig3.pdf", height=7, width=7)
#
ggplot(df5, aes(x=sophistication, y=Percentage, fill=Network)) + ylim(0,100.1) + facet_wrap(~presvote) +
geom_bar(stat="identity") + scale_fill_viridis(discrete=TRUE) + xlab("Political information") +
ggtitle("Composition of voters' discussion networks\n2016 Cooperative Congressional Election Study (weighted)") +
theme(plot.title = element_text(hjust = 0.5))
#
dev.off()
#
# County partisanship
#
df6 <- ddply(df2,.(presvote, countytype),
    function(x) with(x,
      data.frame(100*wtd.table(networkmix, weights=weight)/sum(wtd.table(networkmix, weights=weight)))))
#
names(df6)[names(df6)=="Var1"] <- "Network"
names(df6)[names(df6)=="Freq"] <- "Percentage"
#
df6 <- na.omit(df6)
#
pdf("Dropbox/networks/images/anes2000_fig4.pdf", height=7, width=7)
#
ggplot(df6, aes(x=countytype, y=Percentage, fill=Network)) + ylim(0,100.1) + facet_wrap(~presvote) +
geom_bar(stat="identity") + scale_fill_viridis(discrete=TRUE) + xlab("County type") +
ggtitle("Composition of voters' discussion networks\n2016 Cooperative Congressional Election Study (weighted)") +
theme(plot.title = element_text(hjust = 0.5))
#
dev.off()
#
#
#
#
#    BAYESIAN ALDRICH-MCKELVEY SCALING
#  **** Use only standard 7-pt data ****
#
selfplacement <- V000440
selfplacement[selfplacement < 1 | selfplacement > 7] <- NA
#
placements <- data.frame(V001382, V001383, V000449, V000456, V000466)
placements[placements < 1 | placements > 7] <- NA
names(placements) <- c("DemParty","RepParty","Clinton","Gore","Bush")
#
combined.placements <- data.frame(self=selfplacement, placements)
#
U <- as.matrix(combined.placements)
U[is.na(U)] <- 999
AM.result <- aldmck(U, polarity=2, respondent=1, missing=999, verbose=TRUE)
#
cutoff <- 3   # This section removes respondents who placed less than 3 stimuli, with at least two unique valid placements
index <- rowSums(!is.na(placements)) >= cutoff & apply(placements, 1, function(x){length(na.omit(unique(as.numeric(x))))}) > 1
#
T <- combined.placements[index,]
#
T <- T - 4
#
self <- T[,1]
TT <- T[,-1]
#
N <- nrow(TT)
q <- ncol(TT)
z <- TT
#
modelstring <-
"model{
for(i in 1:N){   ##loop through respondents
        for(j in 1:q){ ##loop through stimuli
			z[i,j] ~ dnorm(mu[i,j], tau[i,j])
			mu[i,j] <- a[i] + b[i]*zhat[j]
			tau[i,j] <- tauj[j] * taui[i]
}
}
##priors on a and b
for(i in 1:N){
a[i] ~ dunif(-100,100)
b[i] ~ dunif(-100,100)
}
##priors on variance
for(j in 1:q){
tauj[j] ~ dgamma(1,.1)
}
for(i in 1:N){
taui[i] ~ dgamma(ga,gb)
}
ga ~ dgamma(1,.1)
gb ~ dgamma(1,.1)
##priors on zhat
zhat[1] ~ dnorm(0,1)T(-1.1,-0.9)
zhat[2] ~ dnorm(0,1)T(0.9,1.1)
for (j in 3:q){
zhat[j] ~ dnorm(0,1)
}
sigmai.sq <- 1/taui
sigmaj.sq <- 1/tauj
}"
#
inits <- function(){
list(
a=runif(N,-2,2),
b=runif(N,0.5,2),
zhat=c(-1,1,rnorm(q-2,0,0.5))
)}
#
#BAM.sim <- jags.model(textConnection(modelstring),
#    data = list('z' = z, 'q' = q, 'N' = N),
#    inits = inits, n.chains = 2, n.adapt = 2500)
#
#update(BAM.sim, n.iter=2500)
#
params <- c("zhat", "a", "b")
#samps <- coda.samples(BAM.sim, params, n.iter=5000, thin=5)
#
# Write results to disk:
#save(samps,file="Dropbox/networks/analysis/ANES2000_samps.Rda")
#
# Load results from disk:
load("Dropbox/networks/analysis/ANES2000_samps.Rda")
#
libcon_a <- samps[,grep(pattern="a[", colnames(samps[[1]]), fixed=TRUE)]
libcon_b <- samps[,grep(pattern="b[", colnames(samps[[1]]), fixed=TRUE)]
libcon_all.a <- do.call("rbind", libcon_a)
libcon_all.b <- do.call("rbind", libcon_b)
#
alpha.mu <- summary(libcon_a)[[1]][,1]
beta.mu <- summary(libcon_b)[[1]][,1]
#
# Note: we need to map back to original data
#
alpha.means <- NA
alpha.means[index] <- alpha.mu
#
beta.means <- NA
beta.means[index] <- beta.mu
#
alpha.means[beta.means < 0] <- NA
#
# Calculate ideal points
#
libcon_self <- self
#
libcon_nsamp <- nrow(libcon_all.a)
libcon_nresp <- ncol(libcon_all.a)
#
libcon_idealpt <- rep(NA, libcon_nsamp*libcon_nresp)
dim(libcon_idealpt) <- c(libcon_nsamp,libcon_nresp)
for (i in 1:libcon_nresp){
for (j in 1:libcon_nsamp){
libcon_idealpt[j,i] <- ((libcon_self[i] - libcon_all.a[j,i]) / libcon_all.b[j,i])
}}
#
libcon_idealpt.percentiles <- rep(NA,libcon_nresp*3)
dim(libcon_idealpt.percentiles) <- c(libcon_nresp,3)
colnames(libcon_idealpt.percentiles) <- c("5%","50%","95%")
for (i in 1:libcon_nresp){
libcon_idealpt.percentiles[i,] <- quantile(libcon_idealpt[,i],probs=c(0.05,0.5,0.95),na.rm=TRUE)
}
libcon_idealpt.estimates <- libcon_idealpt.percentiles[,2]
#
bamidealpoints <- NA
bamidealpoints[index] <- libcon_idealpt.estimates
#
#
#
#
df <- data.frame(alpha.means, bamidealpoints, party=factor(partyid, labels=c("D","I","R")), networkmix, weight)
df2 <- na.omit(df)
#
#dfAtt<-ddply(df2,~party+networkmix,function(x) mean_cl_boot(x$alpha.means, conf.int=.8))
#
#
A1 <- ggplot(df2, aes(x=networkmix, y=alpha.means, fill=networkmix, weight=weight)) +
geom_boxplot() +
#scale_fill_brewer(palette="Spectral") +
scale_fill_viridis(discrete=TRUE, option="viridis") +
xlab("Network Composition") +
ylab(expression(paste("Distortion Parameter (", alpha[i], ")"))) +
ggtitle("Ideological Bias\n") +
guides(fill=guide_legend(title="")) +
theme(plot.title = element_text(hjust = 0.5),
#		axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
theme(legend.position="none")
#theme(plot.title = element_text(hjust = 0.5), legend.position="none")
#
A2 <- ggplot(df2, aes(x=networkmix, y=bamidealpoints, fill=networkmix, weight=weight)) +
geom_boxplot() +
#scale_fill_brewer(palette="Spectral") +
scale_fill_viridis(discrete=TRUE, option="viridis") +
xlab("Network Composition") +
ylab(expression(paste("Corrected Liberal-Conservative Self-Placement (", x[i], ")"))) +
ylim(-5,5) +
ggtitle("Ideological Ideal Point\n") +
guides(fill=guide_legend(title="")) +
theme(plot.title = element_text(hjust = 0.5),
#		axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())
#theme(plot.title = element_text(hjust = 0.5), legend.position="none")
#
pdf("Dropbox/networks/images/ANES2000_BAM_alphaidealpoints_networks.pdf", height=5, width=10)
grid.arrange(A1, A2, ncol=2, widths=c(4.1,5.9))
dev.off()
#
#
#
#
# DISTORTION PARAMETER
#
temp.df <- data.frame(scale(alpha.means), relevel(networkmix, ref = 3),
	rescale01(partyid7), rescale01(ideology), rescale01(income), rescale01(age), rescale01(education),
	rescale01(female), rescale01(black), rescale01(latino), rescale01(unionmember),
	rescale01(religiosity.scores), county.gop2000, politicalsophistication.scores)
colnames(temp.df) <- c("alphameans", "networkmix",
	"partyid", "ideoselfplace", "income", "age", "education",
	"female", "black", "latino", "unionmember",
	"religiosity", "countygop2000", "politicalsophistication")
#
temp.df2 <- na.omit(temp.df)
#
#temp.df2$politicalsophistication <- factor(ntile(temp.df2$politicalsophistication, 2),
#   labels=c("Below Median","Above Median"))
#
temp.df2$politicalsophistication3[temp.df2$politicalsophistication < 0.34] <- 1
temp.df2$politicalsophistication3[temp.df2$politicalsophistication > 0.34 & temp.df2$politicalsophistication < 0.68] <- 2
temp.df2$politicalsophistication3[temp.df2$politicalsophistication > 0.68] <- 3
#
networkmix.binary <- model.matrix(~ (temp.df2$networkmix + temp.df2$politicalsophistication)^2)
temp.df2$onlygore <- networkmix.binary[,2]
temp.df2$gorenotsure <- networkmix.binary[,3]
temp.df2$bushnotsure <- networkmix.binary[,4]
temp.df2$onlybush <- networkmix.binary[,5]
temp.df2$onlygore.sophistication <- networkmix.binary[,7]
temp.df2$gorenotsure.sophistication <- networkmix.binary[,8]
temp.df2$bushnotsure.sophistication <- networkmix.binary[,9]
temp.df2$onlybush.sophistication <- networkmix.binary[,10]
#
temp.df2$partyidstrength <- NA
temp.df2$partyidstrength <- abs(temp.df2$partyid - 0.5)*2
#
onlygore.interflexplot <- inter.binning(Y = "alphameans", D = "onlygore", X = "politicalsophistication", Z = c("partyid", "ideoselfplace",
	"income", "age", "education", "female", "black", "latino", "unionmember", "religiosity", "countygop2000",
	"gorenotsure","bushnotsure","onlybush","gorenotsure.sophistication","bushnotsure.sophistication","onlybush.sophistication"),
	cutoffs=c(0,0.25,0.5,0.75,1), data = temp.df2, Ylabel = "Alpha", Dlabel = "Only gore", Xlabel="Political Knowledge")
#
gorenotsure.interflexplot <- inter.binning(Y = "alphameans", D = "gorenotsure", X = "politicalsophistication", Z = c("partyid", "ideoselfplace",
	"income", "age", "education", "female", "black", "latino", "unionmember", "religiosity", "countygop2000",
	"onlygore","bushnotsure","onlybush","onlygore.sophistication","bushnotsure.sophistication","onlybush.sophistication"),
	cutoffs=c(0,0.25,0.5,0.75,1), data = temp.df2, Ylabel = "Alpha", Dlabel = "gore + Not Sure", Xlabel="Political Knowledge")
#
bushnotsure.interflexplot <- inter.binning(Y = "alphameans", D = "bushnotsure", X = "politicalsophistication", Z = c("partyid", "ideoselfplace",
	"income", "age", "education", "female", "black", "latino", "unionmember", "religiosity", "countygop2000",
	"onlygore","gorenotsure","onlybush","onlygore.sophistication","gorenotsure.sophistication","onlybush.sophistication"),
	cutoffs=c(0,0.25,0.5,0.75,1), data = temp.df2, Ylabel = "Alpha", Dlabel = "bush + Not Sure", Xlabel="Political Knowledge")
#
onlybush.interflexplot <- inter.binning(Y = "alphameans", D = "onlybush", X = "politicalsophistication", Z = c("partyid", "ideoselfplace",
	"income", "age", "education", "female", "black", "latino", "unionmember", "religiosity", "countygop2000",
	"onlygore","gorenotsure","bushnotsure","onlygore.sophistication","gorenotsure.sophistication","bushnotsure.sophistication"),
	cutoffs=c(0,0.25,0.5,0.75,1), data = temp.df2, Ylabel = "Alpha", Dlabel = "Only bush", Xlabel="Political Knowledge")
#
A <- inter.plot(onlygore.interflexplot,
	ylab="Marginal Effect on alpha", xlab="Moderator: Political Knowledge", main="Network: Only D\n")
B <- inter.plot(gorenotsure.interflexplot,
	ylab="Marginal Effect on alpha", xlab="Moderator: Political Knowledge", main="Network: D + Not Sure\n")
C <- inter.plot(bushnotsure.interflexplot,
	ylab="Marginal Effect on alpha", xlab="Moderator: Political Knowledge", main="\nNetwork: R + Not Sure\n")
D <- inter.plot(onlybush.interflexplot,
	ylab="Marginal Effect on alpha", xlab="Moderator: Political Knowledge", main="\nNetwork: Only R\n")
#
pdf("Dropbox/networks/images/interflex_ANES2000.pdf", height=10, width=10)
#
ggarrange(A, B, C, D,  labels = c("","","",""), ncol = 2, nrow = 2)
#
dev.off()
#
#
alpha.lm.all <- lm(alphameans ~ networkmix +
	partyid + ideoselfplace + income + age + education +
	female + black + latino + unionmember +
	religiosity + countygop2000 + politicalsophistication + politicalsophistication:networkmix,
	data=temp.df2)
#
apsrtable(alpha.lm.all)
#
#
# Interaction tests
#
alpha.lm.all.inttest <- lm(alphameans ~ networkmix +
	partyid + ideoselfplace + income + age + education +
	female + black + latino + unionmember +
	religiosity + countygop2000 +
	factor(politicalsophistication3==3) + factor(politicalsophistication3==3):partyid + factor(politicalsophistication3==3):ideoselfplace + factor(politicalsophistication3==3):income + factor(politicalsophistication3==3):age + factor(politicalsophistication3==3):education +
	factor(politicalsophistication3==3):female + factor(politicalsophistication3==3):black + factor(politicalsophistication3==3):latino + factor(politicalsophistication3==3):unionmember +
	factor(politicalsophistication3==3):religiosity + factor(politicalsophistication3==3):countygop2000 + factor(politicalsophistication3==3):networkmix,
	data=temp.df2)
#
summary(alpha.lm.all.inttest)
#
#
alpha.lm.low <- lm(alphameans ~ networkmix +
	partyid + ideoselfplace + income + age + education +
	female + black + latino + unionmember +
	religiosity + countygop2000,
	data=temp.df2[temp.df2$politicalsophistication3<=2,])
#
alpha.lm.high <- lm(alphameans ~ networkmix +
	partyid + ideoselfplace + income + age + education +
	female + black + latino + unionmember +
	religiosity + countygop2000,
	data=temp.df2[temp.df2$politicalsophistication3==3,])
#
apsrtable(alpha.lm.low, alpha.lm.high)
#
#
# Check partisan strength
#
alpha.lm.low2 <- lm(alphameans ~ networkmix +
	partyid + partyidstrength + ideoselfplace + income + age + education +
	female + black + latino + unionmember +
	religiosity + countygop2000,
	data=temp.df2[temp.df2$politicalsophistication3<=2,])
#
alpha.lm.high2 <- lm(alphameans ~ networkmix +
	partyid + partyidstrength + ideoselfplace + income + age + education +
	female + black + latino + unionmember +
	religiosity + countygop2000,
	data=temp.df2[temp.df2$politicalsophistication3==3,])
#
apsrtable(alpha.lm.low2, alpha.lm.high2)
#



#
# Equality of coefficients tests
#
equalitycoefs <- function(b1,sd1,b2,sd2){
    numerator <- b1 - b2
    denominator <- sqrt(((sd1^2) + (sd2^2)))
    zscore <- numerator / denominator
    pvalue <- (2*pnorm(abs(zscore), lower.tail = F))
    out <- c(numerator, denominator, zscore, pvalue)
    names(out) <- c("numerator", "denominator", "zscore", "pvalue")
    return(out)
}
#
equalitycoefs(b1=-0.13, sd1=0.21, b2=-0.11, sd2=0.14)
equalitycoefs(b1=-0.09, sd1=0.26, b2=-0.07, sd2=0.20)
equalitycoefs(b1=0.06, sd1=0.26, b2=-0.33, sd2=0.23)
equalitycoefs(b1=-0.23, sd1=0.20, b2=-0.37, sd2=0.13)
#
equalitycoefs(b1=-0.23, sd1=0.20, b2=-0.13, sd2=0.21)
equalitycoefs(b1=-0.37, sd1=0.13, b2=-0.11, sd2=0.14)
#
equalitycoefs(b1=-0.10, sd1=0.29, b2=-0.26, sd2=0.19)
#
#
# Are weights needed? No correlation = non-informative sampling (Pfeffermann & Sverchkov, 1999)
#
cor(alpha.lm.low$residuals, weight[complete.cases(temp.df)][temp.df2$politicalsophistication3<=2])
cor(alpha.lm.high$residuals, weight[complete.cases(temp.df)][temp.df2$politicalsophistication3==3])
#
apsrtable(alpha.lm.low, alpha.lm.high, stars=1, lev=c(0.1))
apsrtable(alpha.lm.low, alpha.lm.high, stars=1, lev=c(0.05))
apsrtable(alpha.lm.low, alpha.lm.high, stars=1, lev=c(0.01))
#
#
alpha.lm.low.weights <- lm(alphameans ~ networkmix +
	partyid + partyidstrength + ideoselfplace + income + age + education +
	female + black + latino + unionmember +
	religiosity + countygop2000,
	data=temp.df2[temp.df2$politicalsophistication3<=2,],
	weights=weight[complete.cases(temp.df)][temp.df2$politicalsophistication3<=2])
#
alpha.lm.high.weights <- lm(alphameans ~ networkmix +
	partyid + partyidstrength + ideoselfplace + income + age + education +
	female + black + latino + unionmember +
	religiosity + countygop2000,
	data=temp.df2[temp.df2$politicalsophistication3==3,],
	weights=weight[complete.cases(temp.df)][temp.df2$politicalsophistication3==3])
#
#
apsrtable(alpha.lm.low.weights, alpha.lm.high.weights, stars=1, lev=c(0.1))
apsrtable(alpha.lm.low.weights, alpha.lm.high.weights, stars=1, lev=c(0.05))
apsrtable(alpha.lm.low.weights, alpha.lm.high.weights, stars=1, lev=c(0.01))
#
#
alphadf.ANES2000 <- data.frame(alpha.means=temp.df2$alphameans,
	networkmix=networkmix[complete.cases(temp.df)],
	politicalsophistication.binary=temp.df2$politicalsophistication,
	weight=weight[complete.cases(temp.df)],
	partyid7=partyid7[complete.cases(temp.df)])
alphadf.ANES2000$politicalsophistication.binary <- ordered(alphadf.ANES2000$politicalsophistication.binary, levels = c("Low Knowledge", "High Knowledge"))
#
#
#
#
X <- na.omit(data.frame(alpha.means, party=factor(partyid, labels=c("Democrat","Independent","Republican")), networkmix))
#
#plot(allEffects(combined.lm), multiline=TRUE)
#
XX <- aggregate(alpha.means ~ party + networkmix, data=X, mean)
XX2 <- aggregate(alpha.means ~ party + networkmix, data=X, length)
#XX <- aggregate(alpha.means ~ party + networkmix, data=X, function(x){sd(x, na.rm=TRUE)/sqrt(length(x))})
#
pdf("Dropbox/networks/images/ANES2000_BAM_alphameans_networks.pdf", height=5, width=6.5)
#
ggplot(data=XX, aes(x=party,y=alpha.means,group=networkmix,color=networkmix)) +
	geom_line(lwd=1.5) +
	scale_color_viridis(discrete=TRUE, option="viridis") +
	xlab("\nParty Identification") +
	ylab(expression(paste("Distortion Parameter (", alpha[i], ")"))) +
	ylim(-0.8,0.8) +
	ggtitle("Ideological Bias (Group Means)\n") +
	guides(color=guide_legend(title="Network Composition")) +
	theme(plot.title = element_text(hjust = 0.5))
#
dev.off()
#
#
#
#
alphadf.ANES2000 <- data.frame(alpha.means, networkmix, politicalsophistication.binary, weight, partyid7)
alphadf.ANES2000$politicalsophistication.binary <- factor(alphadf.ANES2000$politicalsophistication.binary, labels = c("Low Knowledge", "High Knowledge"))
#
save(alphadf.ANES2000, file="Dropbox/networks/analysis/alphadfANES2000.rda")
#load("Dropbox/networks/analysis/alphadfANES2000.rda")
#
biasplot2000 <- ggplot(na.omit(alphadf.ANES2000), aes(x=networkmix, y=alpha.means, fill=networkmix, weight=weight)) +
geom_boxplot() +
facet_wrap(~politicalsophistication.binary) +
#scale_fill_brewer(palette="Spectral") +
scale_fill_viridis(discrete=TRUE, option="viridis") +
xlab("Network Composition") +
ylab(expression(paste("Distortion Parameter (", alpha[i], ")"))) +
ggtitle("Ideological Bias (ANES 2000)\n") +
guides(fill=guide_legend(title="Network Composition")) +
theme(plot.title = element_text(hjust = 0.5),
#		axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())
#theme(legend.position="none")
#theme(plot.title = element_text(hjust = 0.5), legend.position="none")
#
biasplot2000
#
#
